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Abstract: This paper presents preliminary computational aeroelastic analysis results 

generated in preparation for the first Aeroelastic Prediction Workshop (AePW). These 
results were produced using FUN3D software developed at NASA Langley and are com- 
pared against the experimental data generated during the High REynolds Number Aero- 
Structural Dynamics (HIRENASD) Project. The HIRENASD wind-tunnel model was 
tested in the European Transonic Windtunnel in 2006 by Aachen University's Department 
of Mechanics with funding from the German Research Foundation. The computational 
effort discussed here was performed (1) to obtain a preliminary assessment of the ability 
of the FUN3D code to accurately compute physical quantities experimentally measured 
on the HIRENASD model and (2) to translate the lessons learned from the FUN3D anal- 
ysis of HIRENASD into a set of initial guidelines for the first AePW, which includes test 
cases for the HIRENASD model and its experimental data set. This paper compares the 
computational and experimental results obtained at Mach 0.8 for a Reynolds number of 
7 million based on chord, corresponding to the HIRENASD test conditions No. 132 and 
No. 159. Aerodynamic loads and static aeroelastic displacements are compared at two 
levels of the grid resolution. Harmonic perturbation numerical results are compared with 
the experimental data using the magnitude and phase relationship between pressure coef- 
ficients and displacement. A dynamic aeroelastic numerical calculation is presented at one 
wind-tunnel condition in the form of the time history of the generalized displacements. 
Additional FUN3D validation results are also presented for the AGARD 445.6 wing data 
set. This wing was tested in the Transonic Dynamics Tunnel and is commonly used in 
the preliminary benchmarking of computational aeroelastic software. 

1 INTRODUCTION 

The fundamental technical challenge in computational aeroelasticity is the accurate pre- 
diction of unsteady aerodynamic phenomena and their effect on the aeroelastic response 
of a vehicle. Currently, a benchmarking “standard” for use in validating the accuracy of 
computational aeroelasticity codes does not exist. Many aeroelastic data sets have been 
obtained in wind-tunnel and flight testing; however, none have been globally recognized 
as an ideal data set. There are numerous reasons for this. One is that often such aeroe- 
lastic data sets focus on the aeroelastic phenomena alone (flutter, for example) and do 
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not contain associated information, such as unsteady pressures or structural deflections. 
Other available data sets focus solely on the unsteady pressures and do not address the 
aeroelastic phenomena. Other deficiencies can include omission of relevant data, such as 
flutter frequency, or the acquisition of only qualitative deflection data. In addition to these 
content deficiencies, all of the available data sets present both experimental and computa- 
tional technical challenges. Experimental issues include facility influences, nonlinearities 
beyond those being modeled, and data post processing. From the computational per- 
spective, technical challenges include modeling geometric complexities, coupling between 
the flow and the structure, turbulence modeling, grid issues, and boundary conditions. 
An Aeroelasticity Benchmark Assessment task was initiated at NASA Langley in October 
2009 with the primary objectives of (1) examining the existing potential experimental data 
sets and selecting the one(s) viewed as the most suitable for computational benchmarking 
and (2) performing an initial computational evaluation of these configurations using the 
NASA Langley-developed computational aeroelastic (CAE) software FUN3D [1] as part 
of its code validation process. A successful effort will result in the identification of a focus 
problem for government, industry, and academia to use in demonstrating and comparing 
codes, methodologies, and experimental data to advance the state of the art. Ideally, such 
a focus problem would be but the first of many put forth for this purpose, with a future 
goal being the design, fabrication, and testing of an aeroelastic model recognized by the 
community as a benchmark aeroelastic standard. 

Excellent examples of such a progression and escalation of code validation in the interna- 
tional community are the series of four AIAA Drag Prediction Workshops (DPWs) [2] that 
have been held since 2001 and the AIAA High Lift Prediction Workshop (HiLiftPW) [3] 
that was held in 2010. These workshops had three main objectives. The first was to assess 
the ease and practicality of using state-of-the-art computational methods for aerodynamic 
load prediction. The second was to impartially evaluate the effectiveness of the Navier- 
Stokes solvers, and the final objective was to identify areas for improvement. The structure 
of the DPW and the HiLiftPW provides a template for other computational communities 
seeking similar improvements in accuracy within their own fields. The examination and 
selection of aeroelastic data sets within the Aeroelasticity Benchmark Assessment task 
together with the computational evaluation of these configurations led to initiation of an 
Aeroelastic Prediction Workshop (AePW) series [4] . 

The main focus of this paper will be on comparing the computational aeroelastic results 
generated using FUN3D software against the corresponding experimental data acquired 
during the High REynolds Number Aero-Structural Dynamics (HIRENASD) project [5- 
7], which has also been selected to be part of the first AePW. This computational effort 
was performed to obtain a preliminary assessment of the ability of the FUN3D code to 
accurately compute quantities experimentally measured on the HIRENASD model and 
to begin the work of establishing guidelines for the AePW. A secondary intent of this 
paper is to show additional validation of FUN3D using the AGARD 445.6 wing data set, 
focusing on prediction of the flutter boundary. 

Information relevant to the numerical methods employed will be presented first, including 
those for grid generation and rigid steady, static aeroelastic, forced oscillation and dynamic 
aeroelastic analyses. Details associated with the AGARD 445.6 wing and HIRENASD 
analyses will then be presented. 
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2 NUMERICAL METHOD 


2.1 Grid Generation 


Unstructured tetrahedral grids used in this study were generated using VGRID [8] with 
input prepared using GridTool [9]. For the HIRENASD configuration, the boundary layer 
grid was converted into prism elements using preprocessing options within the FUN3D 
software. For the AGARD 445.6 wing grids, the boundary layer consisted of tetrahedral 
elements. Equation 1 describes the grid point distribution normal to the body, where A n 
is the normal spacing of the n th layer, Ao defines an initial cell height, and the variables 
rq and 72 define the geometric growth rate and the exponential growth rate, respectively. 

A„ = A 0 | 1 +r,(l r r 2 )"-'r' (1) 


The off-surface length scale growth rate is governed by the parameters T i and F 2 according 
to the following equations: 
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and As 0 is the primary grid spacing at the source, As| r | is the primary grid spacing at 
the target point, 77 is the hybrid growth exponent, |r| is the euclidean distance between 
a target point and a source, a is the source strength, and As is the characteristic length 
(average mesh spacing) [8]. The magnitudes of the growth rate parameters used in this 
study are identified in specific sections of this paper describing analyzed configurations. 


2.2 Rigid Steady Flow Analysis 

Solutions to the Reynolds-Averaged Navier-Stokes (RANS) equations were computed us- 
ing the FUN3D flow solver. Turbulence closure was obtained using either the Spalart- 
Allmaras [10] one-equation model or the shear-stress transport (SST) model by Menter [11], 
Flux limitation was accomplished with either an augmented van Leer [12] or Venkatakr- 
ishnan [13] limiter. Inviscid fluxes were computed using either the Low-Diffusion Flux- 
Splitting Scheme (LDFSS) of Edwards [14] or the Roe scheme [15]. Unless stated oth- 
erwise, the solutions in this paper were obtained with an augmented van Leer limiter, 
LDFSS for inviscid fluxes, and the Spalart-Allmaras turbulence model. For the asymp- 
totically steady cases under consideration, time integration was accomplished by an Euler 
implicit backwards difference scheme, with local time stepping to accelerate convergence. 
Most of the cases in this study were run for 5000 iterations to achieve convergence of 
forces and moments to within ±0.5% of the average of their last 1000 iterations. 


2.3 Dynamic Aeroelastic Analysis 

Recently, a dynamic aeroelastic capability was added to the FUN3D solver [16]. For 
structural dynamics analysis, FUN3D is now capable of being loosely coupled with an 


3 



external finite element solver [17], or in the case of the linear structural dynamics used in 
this study, an internal modal structural solver can be utilized. This modal solver is formu- 
lated and implemented in FUN3D in a similar manner to other Langley aeroelastic codes 
(CAP-TSD [18] and CFL3D [19]). For the computations presented here, the structural 
modes were obtained via a normal modes analysis (solution 103) with the Finite Element 
Model (FEM) solver MSC Nastran™ [20]. Modes were interpolated to the surface mesh 
using the method developed by Samareh [21]. 

The dynamic analysis was performed in a three-step process. First, the steady CFD so- 
lution was obtained on the rigid vehicle. Next, a static aeroelastic solution was obtained 
by continuing the CFD analysis in a time accurate mode, allowing the structure to de- 
form. A high value of structural damping was used ( 0.99) so the structure could find 
its equilibrium position with respect to the mean flow before the dynamic response was 
started. Finally, for the dynamic response, the damping was set to the expected value 
( 0.00), and the structure was perturbed in generalized velocity for each of the modes 
included in the structural model. The flow was then solved in the time accurate mode. It 
should be noted that a user-specified modal motion is available in FUN3D. In this study, 
for harmonic perturbation, the modal displacement for mode n was computed as: 

q n = A n sin(uj n t ) (4) 

where A n is amplitude, u> n is frequency, and t is time. 

3 AGARD 445.6 WING ANALYSIS 

The AGARD 445.6 wing was tested in the TDT in 1961 [22], Flutter data from this test 
has been publicly available for over 20 years and has been widely used for preliminary com- 
putational aeroelastic benchmarking. The AGARD wing planform was sidewall-mounted 
and had a quarter-chord sweep angle of 45 degrees, an aspect ratio of 1.65, a taper ratio 
of 0.66, a wing semispan of 2.5 feet, and a wing root chord of 1.833 feet. The wing was 
flutter tested in both air and R-12 heavy gas test mediums at Mach numbers from 0.34 to 
1.14 at zero-degrees angle of attack. Unfortunately, this data set lacks unsteady surface 
pressure measurements necessary for more extensive code validation. 

The FUN3D computations for the AGARD 445.6 wing were performed across the entire 
Mach number range of the experimental data, assuming both inviscid and viscous flows 
with air as the working fluid. The first four modes were used in the aeroelastic analysis. 
Historically, the flutter experimental and numerical results for AGARD 445.6 have been 
presented in the form of “flutter speed index” and “frequency ratio”, which are defined 
as: 


flutter speed index = — , frequency ratio = — (5) 

b UJ a y/n UJ a 

where V is the freestream velocity, b is the wing root semichord, p is the mass ratio, u> is 
the flutter frequency, and oj a is the wing first torsion mode frequency. 

The computational process to calculate flutter speed index and the frequency ratio requires 
multiple steps. First, the rigid steady body solution is obtained. The static aeroelstic 
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analysis is not necessary for the wing with the symmetric airfoil at zero angle-of-attack; 
therefore, next, several dynamic aeroelastic FUN3D computations with different values 
of the dynamic pressure (Q) in the vicinity of the experimental flutter dynamic pressure 
are made. The wing response in the form of the time varying generalized variable is then 
computed and used to calculate the damping ratio from the first bending mode generalized 
variable time history using the logarithmic decrement method. For the stable solution, 
the damping ratio is greater than zero, and for the unstable solution, the damping ratio 
is less than zero. The damping ratio and the dynamic pressure are interpolated, and at 
zero damping ratio the dynamic pressure is considered to be the flutter dynamic pressure. 
Once this flutter dynamic pressure is identified, the corresponding flutter frequency is 
determined via an interpolation of the frequencies at that condition. The computational 
points and zero damping ratio interpolation process for the Mach = 0.9 and angle of 
attack a = 0° condition assuming inviscid flow are shown in Figure 1. 
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Figure 1: Determination of the AGARD 445.6 wing flutter dynamic pressure via interpolation at 
Mach = 0.9 and a = 0° assuming inviscid flow in air. 

Figures 2 (a) and (b) show comparisons among the experimental flutter speed index and 
frequency ratio values, respectively, and those obtained using FUN3D and those published 
in the literature [23-26]. In general, in the subsonic flow regime, the computational data 
matches the experimental data well, while a broad range in the computational data is 
observed in the transonic and supersonic flow regimes. The FUN3D results in Figure 2 
emphasize the importance of the viscous flow assumption near the transonic dip and the 
necessity to use a fine grid in the supersonic flow regime. 

The largest discrepancies between the experimental and the numerical data in Figure 2 
occur at the Mach = 1.14 flow condition. Consequently, a grid refinement study was 
conducted at this condition. Figures 3(a) through 3(d) show the planform and Y = 1 
plane views for the baseline and fine grids used in this analysis. These two tetrahedral 
grids are not in a family of successive grid refinements because of differences in the viscous 
layer growth rates (as defined in Section 2.1) and local grid refinements. For the baseline 
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Figure 2: Experimental and computational flutter speed index and frequency ratio versus Mach number 
for the AGARD 445.6 wing. 


grid, the geometric and exponential growth rates were r± — r 2 = 0.08, while for the fine 
grid they were rq = 0.08 and r 2 = 0.00. The baseline grid consisted of four million points, 
and the fine grid consisted of 8 million points. Both grids used the same initial cell height, 
corresponding to y + = 0.3, and the length-scale growth rate parameters (Ti and r 2 ) were 
set to 0.15. Figures 3(e) through 3(h) show surface C p and Mach contours in the Y = 
1 plane. Overall, the grid refinement resulted in better definition of the supersonic flow 
around the wing, Figure 3(h), and significantly improved the flutter boundary prediction. 
Another grid refinement step will be conducted to further verify these results. 
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(a) Baseline surface grid 


(b) Fine surface grid 



(e) Surface C p for the baseline grid 



(g) Mach contours in the Y = 1 plane 
for the baseline grid 


(f) Surface C v for the fine grid 



(h) Mach contours in the Y = 1 plane 
for the fine grid 


Figure 3: Baseline and fine grids, surface C p , and Mach contours in the Y = 1 plane from FUN3D 
solutions at Mach = 1.14 freestream for the AGARD 445.6 wing. 
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4 HIGH REYNOLDS NUMBER AERO-STRUCTURAL DYNAMICS 
(HIRENASD) PROJECT 

The High REynolds Number Aero-Structural Dynamics (HIRENASD) Project [5-7] was 
led by Aachen University’s Department of Mechanics (LFM) with funding from the Ger- 
man Research Foundation (DFG). It was initiated in 2004 to produce a high-quality tran- 
sonic aeroelastic data set at realistic flight Reynolds numbers for a large transport-type 
wing/body configuration and tested in the European Transonic Windtunnel (ETW) in 
2006. The HIRENASD wing planform, shown in Figure 4, is a ceiling-mounted, semi-span, 
clean- wing configuration with a leading-edge sweep of 34 degrees, a span of approximately 
1.3 meters, and a mean aerodynamic chord of 0.3445 meters. It consists of three sections. 
The two outboard sections use an 11-percent thick BAC3-11/RES/30/21 supercritical 
airfoil. The inboard section uses the same airfoil thickened linearly from 11-percent at 
its outer edge to 15-percent at the root. To minimize boundary layer interference during 
testing, a generic fuselage was included. It extended 0.09 meters from the tunnel ceiling 
and was mechanically isolated from the wing by a labyrinth seal. Extensive measure- 
ments were acquired during testing of the HIRENASD model. Instrumentation included 
a six-component balance, Surface Pattern Tracking (SPT) optical markers for surface de- 
formation measurements on the pressure side of the wing, 11 accelerometers, 28 strain 
gages, and 259 unsteady pressure transducers. The pressure transducers were distributed 
along the upper and lower surfaces at seven span sections. 



Figure 4: HIRENASD wing model planform, assembly, and ETW installation photo. Dimensions shown 
are in millimeters. 


The HIRENASD test matrix consisted of both static and dynamic measurements at dif- 
ferent flow conditions, with variations of Reynolds number from 7 million up to 73 million 
based on the mean aerodynamic chord and dynamic pressures up to 130,000 pascals at 
six transonic Mach numbers: 0.70, 0.75, 0.80, 0.83, 0.85, 0.88. The test medium was 


nitrogen. For static testing, pressure distribution and lift and drag polars were acquired 
at various angles of attack. Dynamic testing involved forced vibrations of the wing at the 
natural frequencies of the first bending, second bending, and first torsion modes and was 
performed over the range of Reynolds numbers at different angles of attack. 

The focus of this paper is on comparing the computational results obtained using FUN3D 
with the HIRENASD experimental test cases No. 132 for angle of attack sweep and 
No. 159, where the model was harmonically excited at the second bending mode frequency. 
The rest of this section is organized in the following way. First, the description of the 
grids used in the analysis is presented. Next, the finite element model and the issues 
encountered with interpolating the mode shapes into the surface mesh are described. 
The rigid body and the static aeroelastic results are then compared together against the 
experimental data. Next, the harmonic excitation experimental data are compared against 
the unsteady computational results. Finally, some computational dynamic aeroelastic 
results are briefly introduced. 

4.1 Unstructured Grids 

Three unstructured tetrahedral grids with prism elements in the boundary layer were 
used in this study. The grid parameters, which are based on a Reynolds number of 23.5 
million and c re f = 0.3445 meter, are shown in Table 1. Grids A and B belong to the 
same family of grids. Grid C was constructed to investigate the effect of the wind-tunnel 
ceiling modeled as a viscous surface. In this grid, the location of the outer boundary 
in the computational domain was reduced from 100 to 30 times the c re f. This change 
was arbitrary since no information regarding boundary layer growth inside the ETW was 
readily available. The surface grid extracted from grid A is shown in Figure 5. 



Grid A 

Grid B 

Grid C 

Chordwise Spacing at Wing LE 

0.20%cy oo t 

0.07%c root 

0.07%c roo t 

Chordwise Spacing at Wing TE 

O.lOVoCroot 

0.07%c root 

0.07%c root 

Average Cell y + 

0.444 

0.296 

0.296 

Prism Layer Cells 

30 

30 

30 

Viscous Wall Spacing Ao (meter) 

1.96FT 7 

1.31A- 7 

1.31FT 7 

Grid Size (nodes in millions) 

10 

30 

30 

Outer Boundary 

100c re / 

100c r e / 

30c r ef 

Viscous Stretching r \^2 

0.02, 0.15 

0.02, 0.15 

0.02, 0.15 

ri,r 2 

0.15, 0.15 

0.15, 0.15 

0.15, 0.15 

Wind-Tunnel Ceiling Boundary Cond. 

Symmetry 

Symmetry 

Viscous 


Table 1: VGRID parameters for HIRENASD configuration grids. 


4.2 Finite Element Model 

Two different Finite Element Models (FEMs) are available from the HIRENASD web- 
site [27]. Both are modeled with uniform solid elements. One FEM uses NASTRAN 
hexagonal elements with over 200,000 grid points, while the other uses NASTRAN tetra- 
hedral elements with approximately 170,000 grid points. The two FEMs yield slightly 
different modal frequencies. However, these differences are small (less than 1.3 percent), 
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Figure 5: Grid A surface grid used in FUN3D Figure 6: FEM with hexagonal elements used in 

analysis for the HIRENASD configura- modal analysis for the HIRENASD config- 

tion. uration. 


and the first ten modes are nearly identical. Modal analyses using the available FEMs 
were performed, where the material properties were varied to investigate their impact on 
modal frequencies. These frequencies were plotted (not shown here) against dynamic pres- 
sure normalized by Young’s Modulus, which was a function of temperature. The resulting 
variations of the modal frequencies and mode shapes were assessed to be negligible. For 
the aeroelastic analyses presented in this paper, the first ten modes extracted from the 
NASTRAN model with hexagonal elements, shown in Figure 6, were used. 

4.3 Rigid Body and Static Aeroelastic Analyses 

The rigid body solutions corresponding to HIRENASD test case No. 132 at Mach = 0.8 
and Re = 7 million were obtained at five angles of attack for grids A and B (-1.5, 0.0, 1.5, 
3.0, and 4.5 deg) and one angle of attack for grid C (1.5 deg). The rigid body solutions 
were then used as initial conditions for the corresponding static aeroelastic solutions, 
which in turn were obtained by running FUN3D in the unsteady flow mode with the 
linear structural solver, where the structural damping was set to a large value. The 
solution process required that the mode shapes be interpolated into the surface grid. In 
the first attempt of the static aeroelastic analysis, the mode shapes were interpolated into 
a surface grid consisting of the wing geometry only, even though the wing was attached to 
the fuselage. In this process, FUN3D failed in the beginning of the unsteady flow solution 
due to the negative volumes at the wing-fuselage junction. It was deduced that this 
failure was caused by the discontinuity in mode shape at this location, resulting from the 
geometric differences between structural model and the wing-only shape. A new method 
was subsequently developed, where initially the modes shapes were interpolated into both 
the wing and fuselage. The mode shapes were then set to zero on the fuselage, with the 
exception of a very narrow region near the wing and fuselage junction. In this region, the 
mode shapes were linearly interpolated from zero to the nominal value at the root of the 
wing. This second method resulted in successful execution of the flow solver. Figure 7 
depicts an example of the mode shape interpolation process from FEM to CFD surface 
mesh for the second bending mode (mode 2). Figure 7(a) shows mode 2 as extracted 
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from the FEM, Figure 7(b) shows the interpolation of mode 2 onto the wing-only surface, 
which resulted in flow solver failure due to the negative volumes, and Figure 7(c) shows 
the interpolation of mode 2 onto the wing and fuselage, which resulted in the successful 
execution of FUN3D. 



(a) Mode 2 from FEM (b) Mode 2 interpolation to wing-only, 

which resulted in negative volumes 



(c) Mode 2 interpolation to wing and 
fuselage, resulting in flow solver success- 
ful execution 


Figure 7: Example of the mode shape interpolation process from FEM into CFD surface mesh for the 
HIRENASD configuration. 


The aerodynamic coefficients obtained from rigid body and static aeroelastic calculations 
for grids A, B, and C are shown in Table 2. A decrease in computed drag coefficient is 
observed between grid A and grids B and C. Also, a small lift decrement exists between 
the rigid and deformed (static aeroelastic) configurations. The wing tip displacements 
calculated using grids A and B compared with the experimental SPT data are shown in 
Figure 8. The experimental and computational data compare very well. The wing twist 
experimental SPT data was not available when this document was written; therefore, a 
corresponding plot of the experimental versus computational spanwise wing twist is not 
presented. However, the computational wing twist values (Aa) at the wing tip are shown 
in Table 3. 

As previously mentioned, the HIRENASD experimental unsteady pressure data was col- 
lected at seven span sections, which are identified in Figure 9. Figure 10 presents the C p 
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Grid A 

Steady Rigid 

Solution 

Grid A 

Static Aeroelastic 

Solution 

AoA (°) 

C L 

Cd 

Cm 

C L 

Cd 

Cm 

-1.5 

0.0145 

0.0137 

-0.0964 

0.0099 

0.0138 

-0.0887 

0.0 

0.1811 

0.0132 

-0.3178 

0.1703 

0.0131 

-0.2997 

1.5 

0.3542 

0.0173 

-0.5516 

0.3373 

0.0166 

-0.5231 

3.0 

0.5218 

0.0260 

-0.7770 

0.4996 

0.0242 

-0.7399 

4.5 

0.6493 

0.0408 

-0.9266 

0.6373 

0.0372 

-0.9118 


Grid B 

Steady Rigid 

Solution 

Grid B 

Static Aeroelastic 

Solution 

AoA (°) 

C L 

C D 

Cm 

C L 

C D 

C M 

-1.5 

0.0135 

0.0117 

-0.0974 

0.0087 

0.0117 

-0.0894 

0.0 

0.1806 

0.0113 

-0.3204 

0.1693 

0.0111 

-0.3012 

1.5 

0.3544 

0.0154 

-0.5574 

0.3366 

0.0147 

-0.5271 

3.0 

0.5190 

0.0237 

-0.7794 

0.4964 

0.0219 

-0.7416 

4.5 

0.6304 

0.0373 

-0.8954 

0.6236 

0.0337 

-0.8927 


Grid C 

Steady Rigid 

Solution 

Grid C 

Static Aeroelastic 

Solution 

AoA (°) 

C L 

C D 

Cm 

C L 

C D 

C M 

1.5 

0.3469 

0.0146 

-0.5461 

0.3294 

0.0140 

-0.5163 


Table 2: Computed aerodynamic coefficients for grids A, B, and C (rigid and deformed solutions): 
HIRENASD test case No. 132, Mach = 0.8, Re = 7 million. 



Accelerometer 15 "section 7 ~ 
Section 6 


Section 5 
Section 4 - 


Section 3 
Section 2 — 

Section 1 


Figure 8: Wing tip displacements for HIRENASD Figure 9: HIRENASD span sections 1-7 and 
test case No. 132, Mach = 0.8, accelerometer 15. 

Re = 7 million. 


results at sections 1, 2, 5, and 7 for grids A, B, and C obtained from the rigid steady 
and static aeroelastic solutions. As shown in Figures 10(a) and (b) at section 1 for grids 
A and B, the computed C p values have a very similar distribution to the experimental 
data, but the computed shock appears to be stronger and further aft on the wing for both 
grids. The grid C solution shown in Figure 10(c) captures the shock location very well. 
This suggests that modeling the tunnel wall as a viscous surface and thus capturing the 
boundary layer does have some influence on the wing surface pressure near the ceiling. 
The computed C p at section 2 matches the experimental data quite well for all three 
grids, as shown in Figures 10(d)- (f). The largest discrepancy between the experimental 
and computational data is observed at section 5. Neumann [28] showed that the shape of 
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Grid A 

Static Aero. 

Solution 

Grid B 

Static Aero. 

Solution 

AoA (°) 

A z LE (mm) 

A zte (mm) 

A« (°) 

A z LE (mm) 

A zte (mm) 

Aa (°) 

-1.5 

1.12 

1.93 

-0.32 

1.51 

1.90 

-0.15 

0.0 

7.21 

7.64 

-0.17 

6.92 

7.71 

-0.30 

1.5 

12.89 

13.74 

-0.33 

12.78 

14.02 

-0.48 

3.0 

18.43 

19.68 

-0.49 

18.42 

20.07 

-0.64 

4.5 

22.43 

23.93 

-0.59 

21.69 

23.51 

-0.71 


Table 3: Computed leading- and trailing-edge displacements and the wing twist for grids A and B: HIRE- 
NASD test case No. 132, Mach =0.8, Re = 7 million. 


the trailing edge of the HIRENASD wing has a significant effect on C p at section 5. He 
compared computational results using the Tau [29] code where the trailing edge of the 
wing was considered to be either sharp, blunt, or blunt with rounded corners, which was 
the shape used to generate grids for the FUN3D computations. Neumann’s results for the 
blunt trailing edge with rounded corners (not shown here) match the solutions presented 
in Figures 10(g)-(i). 

4.4 Harmonic Excitation Analysis 

The HIRENASD experimental harmonic excitation tests were conducted to measure the 
interaction between aerodynamics and modal excitations. For test case No. 159 (Mach = 
0.8, Re = 7 million), the model was excited at 78.9 Hz, which closely corresponds to the 
second bending mode frequency. Numerically, the modal excitation is accomplished using 
Equation 4. The unsteady solutions were restarted from the static aeroelastic solutions 
and were run for six cycles, with 64 time steps per cycle and 25 subiterations per time 
step. A drop in the residuals of three to five orders of magnitude was obtained within 
25 subiterations. The unsteady surface pressure and surface shape were collected at 
each time step for four cycle. The results were post-processed to produce the transfer 
function estimate of pressure due to displacement at the specified locations. For this 
study, the location of the displacement used in the analysis coincided with the location of 
accelerometer 15, which is shown in Figure 9. The displacement value was normalized by 
the reference chord. Figures 11 and 12 show the resulting C p magnitude and phase plots, 
respectively, for sections 1 through 7 for grids A, B, and C. For these calculations, the 
peak-to-peak amplitude at accelerometer location 15 was 8.5 mm. The shock locations at 
sections 1-5 are predicted quite well, but the prediction is poor near the wing tip. Figure 13 
shows the C p magnitude plots where the computational solutions were obtained with grid 
A for three different excitation amplitudes: Ampl = 4.25mm, Amp2 = 8.5mm, and Amp3 
= 85mm. Clearly, the computational excitation amplitude has an effect on the numerical 
results. Similarly, Figure 14 shows the magnitude plots where three different methods 
within FUN3D were used to obtain the solution on grid A: (1) LDFSS for inviscid fluxes 
with augmented van Leer limiter and Spalart-Allmaras turbulence model (labeled as ’Grid 
A’), (2) LDFSS for inviscid fluxes with augmented van Leer limiter and SST turbulence 
model (labeled as ’Grid A SST’), and (3) Roe scheme with Venkatakrishnan limiter and 
Spalart-Allmaras turbulence model (labeled as ’Grid A Roe SA’). The results show that 
the solution is more dependent on the choice of the inviscid flux method or limiter than it 
is on the choice of turbulence model. Future analysis will address the effects of the time 
step size on the surface pressure calculation. 
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(a) Grid A Section 1 (b) Grid B Section 1 (c) Grid C Section 1 





(d) Grid A Section 2 (e) Grid B Section 2 (f) Grid C Section 2 





(g) Grid A Section 5 (h) Grid B Section 5 


(i) Grid C Section 5 





(k) Grid B Section 7 


(1) Grid C Section 7 


Figure 10: 


Rigid steady and static aeroelastic surface pressure comparison for grids A, B, and C at 
sections 1, 2, 5, and 7 for HIRENASD test case No. 132, Mach = 0.8, Re = 7 million. 
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Figure 11: Magnitudes of the transfer function of the pressure coefficient due to the displacement at 
accelerometer 15 for grids A, B, and C; HIRENASD test case No. 159, Mach = 0.8, Re = 7 
million. 
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(e) Section 5 


(f) Section 6 



Figure 12: Phase plots of the transfer function of the pressure coefficient due to the displacement at 
accelerometer 15 for grids A, B, and C; HIRENASD test case No. 159, Mach = 0.8, Re = 7 
million. 
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Figure 13: Magnitudes of the transfer function of the pressure coefficient due to the displacement at 
accelerometer 15 for grid A at three excitation amplitudes; HIRENASD test case No. 159, 
Mach = 0.8, Re = 7 million. 
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Figure 14: Magnitudes of the transfer function of the pressure coefficient due to the displacement at 
accelerometer 15 for grid A for different flow solver options; HIRENASD test case No. 159, 
Mach = 0.8, Re = 7 million. 
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4.5 Dynamic Aeroelastic Analysis 


The dynamic aeroelastic analysis capability in FUN3D was exercised using grid A to 
investigate the dynamic aeroelastic behavior of the HIRENASD configuration at flow 
conditions corresponding to test No. 159. In this simulation, the dynamic analysis was re- 
started from the static aeroelastic solution by setting the structural damping to zero value 
and perturbing the structure in generalized velocity for each of the modes used. Ten modes 
were perturbed with the same value of the initial generalized velocity. The time step was 
set sufficiently small such that the cycles of the first and tenth modal frequencies were 
resolved into 1205 time steps and 51 time steps, respectively. The generalized coordinates 
for modes 1 through 6 versus time are shown in Figure 15. The solution shows that modes 
1, 2, 4, 5, and 6 are damped, while mode 3 is neutrally damped. 




•pyu v 


0.2 0.3 

time (seconds) 


Mode 1 
Mode 2 
Mode 3 
Mode 4 
Mode 5 
Mode 6 


0.4 


0.5 


Figure 15: Generalized coordinate as a function of time for modes 1-6, HIRENASD configuration, 
Mach = 0.8, Re = 7 million. 


5 CONCLUDING REMARKS 

Reynolds- Averaged Navier-Stokes analyses, including static aeroelastic, forced excitation, 
and dynamic aeroelastic computations, were performed as part of the validation process 
of the NASA Langley FUN3D software for two configurations: the AGARD 445.6 wing 
and HIRENASD. The AGARD 445.6 wing analysis showed that FUN3D could predict 
the flutter boundary well for this planform. The analysis also showed the importance of 
a viscous flow assumption to accurately predict the flutter boundary at transonic speeds 
near the transonic dip. In addition, a highly resolved grid was required to match the 
experimental data at the supersonic speeds. 

The HIRENASD configuration analysis showed that FUN3D predictions matched the ex- 
perimental data well. The computed static aeroelastic wing tip displacement and the 
magnitude and phase of the transfer function of the pressure coefficient due to normalized 
deflection showed very satisfactory results. Additional grid convergence and computa- 
tional time step sensitivity iterations are needed. 

Lessons learned from the analysis of the HIRENASD configuration have resulted in several 
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recommendations for the first Aeroelastic Prediction Workshop, with the goal of reducing 
uncertainty across participant-generated results. First, common gridding guidelines that 
establish the definition of a grid family need to be developed. A common finite element 
model and a common method of interpolating modes to the surface grid must also be 
established. To avoid ambiguities in the surface definition, such as the wing’s trailing- 
edge geometry, a common baseline geometry is needed. And finally, a common method 
for post-processing the time-accurate data needs to be defined. 
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